An analysis of amaranth data

Visualizing the revised data

The plots below each visualize the prepared data, with colors indicating either accession, species, color, or origin, as well as linear trends per accession. Hover your pointer over the data points for more information about them. Click legend elements once to turn them on or off; double-click to isolate an element; double-click again to reveal hidden elements.

(
  am_data %>%
    dplyr::mutate(`Accession` = stringr::str_c(`Accession ID`,": ",Accession),
                  `Accession` = factor(`Accession`, levels = unique(`Accession`))) %>%
    ggplot(mapping = aes(x = `Diameter (mm)`,
                         y = `Thickness (mm)`,
                         group = `Accession`,
                         color = `Accession`,
                         text = stringr::str_c("Accession: ",`Accession`,"\n",
                                               "Origin: ",`Origin`,"\n",
                                               "Species: ",`Species`,"\n",
                                               "Color: ",`Color`,"\n",
                                               "Diameter: ",`Diameter (mm)`," mm \n",
                                               "Thickness: ",`Thickness (mm)`," mm"))) +
    geom_smooth(method=lm) +
    geom_point()
) %>%
  plotly::ggplotly(tooltip = "text")

(
  am_data %>%
    ggplot(mapping = aes(x = `Diameter (mm)`,
                         y = `Thickness (mm)`,
                         group = `Accession ID`,
                         color = `Species`,
                         text = stringr::str_c("Accession: ",`Accession`,"\n",
                                               "Origin: ",`Origin`,"\n",
                                               "Species: ",`Species`,"\n",
                                               "Color: ",`Color`,"\n",
                                               "Diameter: ",`Diameter (mm)`," mm \n",
                                               "Thickness: ",`Thickness (mm)`," mm"))) +
    geom_smooth(method=lm) +
    geom_point()
) %>%
  plotly::ggplotly(tooltip = "text")

(
  am_data %>%
    ggplot(mapping = aes(x = `Diameter (mm)`,
                         y = `Thickness (mm)`,
                         group = `Accession ID`,
                         color = `Color`,
                         text = stringr::str_c("Accession: ",`Accession`,"\n",
                                               "Origin: ",`Origin`,"\n",
                                               "Species: ",`Species`,"\n",
                                               "Color: ",`Color`,"\n",
                                               "Diameter: ",`Diameter (mm)`," mm \n",
                                               "Thickness: ",`Thickness (mm)`," mm"))) +
    geom_smooth(method=lm) +
    geom_point()
) %>%
  plotly::ggplotly(tooltip = "text")

(
  am_data %>%
    ggplot(mapping = aes(x = `Diameter (mm)`,
                         y = `Thickness (mm)`,
                         group = `Accession ID`,
                         color = `Origin`,
                         text = stringr::str_c("Accession: ",`Accession`,"\n",
                                               "Origin: ",`Origin`,"\n",
                                               "Species: ",`Species`,"\n",
                                               "Color: ",`Color`,"\n",
                                               "Diameter: ",`Diameter (mm)`," mm \n",
                                               "Thickness: ",`Thickness (mm)`," mm"))) +
    geom_smooth(method=lm) +
    geom_point()
) %>%
  plotly::ggplotly(tooltip = "text")

Pairwise Hotelling’s \(T^2\)

combinations <- expand.grid(am_data$`Accession ID` %>% 
                              unique(),
                            am_data$`Accession ID` %>% 
                              unique()) %>%
  as_tibble()

combinations_char <- combinations %>%
  dplyr::mutate_all(.funs = funs(as.character))

pairwise_hotelling <- combinations_char %>%
  purrrlyr::by_row(..f = function(x){
    test <- Hotelling::hotelling.test(`Diameter (mm)` + `Thickness (mm)`~`Accession ID`,
                                      data = am_data,
                                      pair = c(x$Var1,x$Var2))
    return(list(`$T^2$` = test$stats$statistic,
                m = test$stats$m,
                # df = test$stats$df,
                # nx = test$stats$nx,
                # nx = test$stats$ny,
                # p = test$stats$p,
                p = test$pval))
  },
  .collate = "list") %$%
  .out %>%
  purrr::transpose() %>%
  as.tibble() %>%
  tidyr::unnest() %>%
  cbind(combinations, .)

pairwise_hotelling.char <- pairwise_hotelling %>%
  dplyr::mutate_at(.vars = vars(p), .funs = ~formatC(.,format = "e", digits = 2)) %>%
  dplyr::mutate_at(.vars = vars(`$T^2$`,`m`), .funs = ~sprintf("%7.3f", .))

These data were not collected using a balanced sampling strategy, as would be the case in, say, formal crop trials. Instead, data are esoteric—the contemporary data are drawn from a diverse set of geographic locations under unknown cultivation methods. Furthermore, the contemporary populations analyzed are the result of an unknown number of grow-outs for seed preservation by Native Seeds/SEARCH, which may have affected population structure since accessioning. The ancient seeds, of course, have unknown provenience (in the sense that they could have been transported to the Dyck cliff dwelling from an unknown distance). Because of these restrictions, our dataset is not conducive to a multivariate factorial MANOVA design. The most robust analysis we can hope for is to assess pairwise two-sample tests between the accessions; here, we do so not only between the ancient and modern accessions, but among the modern accessions as well. We use the pairwise significance to assess whether the ancient seeds can be distinguished from the modern accessions, and the test statistic as a distance measure between accessions in a cluster analysis among the accessions.

In the pairwise Hotelling’s \(T^2\) test result tables below, the significant results (\(p < 0.05\)) are denoted in red. Pairs that are displayed in black are not significantly different from one another. Most accessions are statistically distinguishable from one another. The Dyck cliff dwelling cannot be statistically distinguished from accession U (C02-014 Marbled (C016); Morelos, Mexico; A. cruentus; black) along only the diameter and thickness dimensions.

pairwise_hotelling.char %>%
  dplyr::mutate_at(.vars = vars(`$T^2$`:`p`), .funs = ~cell_spec(.,
                                                                 "html",
                                                                 color = ifelse(as.numeric(pairwise_hotelling$p) < 0.05,
                                                                                "red",
                                                                                "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `$T^2$`) %>%
  tidyr::spread(Var1,`$T^2$`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r")

pairwise_hotelling.char %>%
  dplyr::mutate_at(.vars = vars(`$T^2$`:`p`), .funs = ~cell_spec(.,
                                                                 "html",
                                                                 color = ifelse(as.numeric(pairwise_hotelling$p) < 0.05,
                                                                                "red",
                                                                                "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `m`) %>%
  tidyr::spread(Var1,`m`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r")

pairwise_hotelling.char %>%
  dplyr::mutate_at(.vars = vars(`$T^2$`:`p`), .funs = ~cell_spec(.,
                                                                 "html",
                                                                 color = ifelse(as.numeric(pairwise_hotelling$p) < 0.05,
                                                                                "red",
                                                                                "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `p`) %>%
  tidyr::spread(Var1,`p`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r")

Pairwise two-sample kernel density based test

pairwise_kde <- combinations_char %>%
  purrrlyr::by_row(..f = function(x){
    
    d1 <- am_data %>%
      dplyr::filter(`Accession ID` == x$Var1) %>%
      dplyr::select(`Diameter (mm)`,`Thickness (mm)`) %>%
      as.matrix()
    
    d2 <- am_data %>%
      dplyr::filter(`Accession ID` == x$Var2) %>%
      dplyr::select(`Diameter (mm)`,`Thickness (mm)`) %>%
      as.matrix()
    
    
    
    
    test <- ks::kde.test(x1 = d1,
                         x2 = d2)
    
    return(list(T = test$Tstat,
                z = test$zstat,
                p = test$pvalue))
  },
  .collate = "list") %$%
  .out %>%
  purrr::transpose() %>%
  as.tibble() %>%
  tidyr::unnest() %>%
  cbind(combinations, .)

pairwise_kde.char <- pairwise_kde %>%
  dplyr::mutate_at(.vars = vars(p), .funs = ~formatC(.,format = "e", digits = 2)) %>%
  dplyr::mutate_at(.vars = vars(`T`,`z`), .funs = ~sprintf("%7.3f", .))

Kernel density based two-sample tests can provide a more robust comparison between multivariate datasets, as they are nonparametric and asymptotically normal—they compare empirical distributions, as opposed to merely comparing sample means and standard deviations as in parametric tests like Hotelling’s \(T^2\). Here, we use the method developed by Duong, Goud, and Schauer (2012), as provided by the kde.test function in the ks package in R.

The results of the pairwise two-sample kernel density based test are very similar to Hotelling’s \(T^2\); the Dyck samples are distinct from all modern accessions except accession U.

pairwise_kde.char %>%
  dplyr::mutate_at(.vars = vars(`T`:`p`), .funs = ~cell_spec(.,
                                                             "html",
                                                             color = ifelse(as.numeric(pairwise_kde$p) < 0.05,
                                                                            "red",
                                                                            "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `T`) %>%
  tidyr::spread(Var1,`T`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r")

pairwise_kde.char %>%
  dplyr::mutate_at(.vars = vars(`T`:`p`), .funs = ~cell_spec(.,
                                                             "html",
                                                             color = ifelse(as.numeric(pairwise_kde$p) < 0.05,
                                                                            "red",
                                                                            "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `p`) %>%
  tidyr::spread(Var1,p) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r")

Pairwise slope analysis

pairwise_slope <- combinations_char %>%
  purrrlyr::by_row(..f = function(x){
    
    d1 <- am_data %>%
      dplyr::filter(`Accession ID` == x$Var1) %>%
      dplyr::mutate(`Accession ID` = "x")
    
    d2 <- am_data %>%
      dplyr::filter(`Accession ID` == x$Var2) %>%
      dplyr::mutate(`Accession ID` = "y")
    
    m1 <- lm(`Thickness (mm)` ~ `Diameter (mm)` * `Accession ID`,
             data = bind_rows(d1,d2))
    
    m2 <- lm(`Thickness (mm)` ~ `Diameter (mm)` + `Accession ID`,
             data = bind_rows(d1,d2))
    
    test <- anova(m1,m2)
    
    return(list(`Δm` = (lm(formula = `Thickness (mm)` ~ `Diameter (mm)`,
                           data = d1)$coefficients[[2]]) - 
                  (lm(formula = `Thickness (mm)` ~ `Diameter (mm)`,
                      data = d2)$coefficients[[2]]),
                `F` = test$F[[2]],
                p = test$`Pr(>F)`[[2]]))
  },
  .collate = "list") %$%
  .out %>%
  purrr::transpose() %>%
  as.tibble() %>%
  tidyr::unnest() %>%
  cbind(combinations, .)

pairwise_slope.char <- pairwise_slope %>%
  dplyr::mutate_at(.vars = vars(p), .funs = ~formatC(.,format = "e", digits = 2)) %>%
  dplyr::mutate_at(.vars = vars(`Δm`,`F`), .funs = ~sprintf("%7.3f", .))

Our final analysis compares the slope of the linear relationship of thickness and diameter for each accession. We ask the question of whether this relationship may be combined with the raw measurements to further define consistent groups. Here, we do so by testing for the interaction between the one of the dependent variables and the factor (accession) using an analysis of covariance (ANCOVA)—we model seed thickness as a response to seed diameter, with and without interaction with accession. If the model with interaction provides a significantly better fit to the data, the slopes within the groups are significantly different.

As above, significant slope differences (\(Δm\)) are shown in red. The results contrast with the analyses above. The relationship of diameter to thickness in Dyck seeds are most similar to those of accessions A (C02-007 Hopi Red Dye (C002); Hopi, Arizona; A. cruentus; black) and E (C02-013 Alegria (C008); Morelos, Mexico; A. cruentus; white), with larger but still statistically insignificant differences with accessions F, K, and Q. In general, accessions from simliar origins have similar slopes—for example, the slopes all the accessions from Morelos, Mexico are mutually indistinguishable. This suggests a relationship between local environment and seed morphology. Future research to assess the consistency of this pattern should be undertaken.

pairwise_slope.char %>%
  dplyr::mutate_at(.vars = vars(`Δm`:`p`), .funs = ~cell_spec(.,
                                                             "html",
                                                             color = ifelse(as.numeric(pairwise_slope$p) < 0.05,
                                                                            "red",
                                                                            "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `Δm`) %>%
  tidyr::spread(Var1,`Δm`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r")

pairwise_slope.char %>%
  dplyr::mutate_at(.vars = vars(`Δm`:`p`), .funs = ~cell_spec(.,
                                                             "html",
                                                             color = ifelse(as.numeric(pairwise_slope$p) < 0.05,
                                                                            "red",
                                                                            "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `F`) %>%
  tidyr::spread(Var1,`F`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r") 

pairwise_slope.char %>%
  dplyr::mutate_at(.vars = vars(`Δm`:`p`), .funs = ~cell_spec(.,
                                                             "html",
                                                             color = ifelse(as.numeric(pairwise_slope$p) < 0.05,
                                                                            "red",
                                                                            "black"))) %>%
  dplyr::select(Var1,
                Var2,
                `p`) %>%
  tidyr::spread(Var1,p) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  knitr::kable(escape = FALSE,
               row.names = TRUE,
               align = "r")

Clustering by Hotelling’s \(T^2\)

fit <- pairwise_hotelling %>%
  dplyr::mutate(`$T^2$` = `$T^2$` * m) %>%
  dplyr::select(Var1,
                Var2,
                `$T^2$`) %>%
  tidyr::spread(Var1,`$T^2$`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  magrittr::divide_by(.,max(.)) %>%
  as.dist() %>%
  hclust(method="ward.D2")
fit$height %<>%
  magrittr::divide_by(.,max(.))

fit %<>%
  dendro_data(type="rectangle")

fit$labels %<>%
  dplyr::left_join(am_data %>%
                     dplyr::select(`Accession ID`,
                                   `Accession`,
                                   Origin,
                                   Species,
                                   Color) %>%
                     dplyr::distinct(),
                   by = c("label" = "Accession ID")) %>%
  dplyr::mutate(`Accession` = stringr::str_c(`label`,": ",Accession),
                  `Accession` = factor(`Accession`, levels = unique(`Accession`)))

dendro_plot <- ggplot() + 
  geom_segment(data = fit$segments,
               aes(x = x,
                   y = y,
                   xend = xend,
                   yend = yend)) +
  coord_flip() +
  scale_y_reverse(limits = c(1.0,-0.6),
                  breaks = seq(0,1,0.2)) +
  ylab("Scaled Distances")+
  ggplot2::theme(panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        axis.title.x = element_text(colour = NA), 
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
        )

The plots below each visualize the data clustered by Hotelling’s \(T^2\), with label colors indicating either species, color, or origin. The distance measure (Hotelling’s \(T^2\)) is standardized to be between zero and one.

Clusters are strongly oriented by color, reinforcing the relationship between color and seed thickness among domesticated amaranths as suggested above. Apparent slight clustering by species is likely due to the fact that all A. hypochondriacus accessions are white. There is no apparent relationship between place of origin and Hotelling’s \(T^2\) distances. The Dyck amaranth seeds cluster with accession U (C02-014 Marbled (C016)).

dendro_plot +
  geom_text(data = fit$labels,
            aes(x = x,
                y = y - 0.025,
                label = Accession,
                hjust = "outward",
                color = Species),
            size = 3) 

dendro_plot +
  geom_text(data = fit$labels,
            aes(x = x,
                y = y - 0.025,
                label = Accession,
                hjust = "outward",
                color = Color),
            size = 3) 

dendro_plot +
  geom_text(data = fit$labels,
            aes(x = x,
                y = y - 0.025,
                label = Accession,
                hjust = "outward",
                color = Origin),
            size = 3) 

Clustering by \(\lvertΔm\rvert\) — the absolute difference in slopes

fit <- pairwise_slope %>%
  dplyr::mutate(`Δm` = abs(`Δm`)) %>%
  dplyr::select(Var1,
                Var2,
                `Δm`) %>%
  tidyr::spread(Var1,`Δm`) %>%
  magrittr::set_rownames(.,.$Var2) %>%
  dplyr::select(-Var2) %>%
  magrittr::set_rownames(.,colnames(.)) %>%
  magrittr::divide_by(.,max(.)) %>%
  as.dist() %>%
  hclust(method="ward.D2")

fit$height %<>%
  magrittr::divide_by(.,max(.))

fit %<>%
  dendro_data(type="rectangle")

fit$labels %<>%
  dplyr::left_join(am_data %>%
                     dplyr::select(`Accession ID`,
                                   `Accession`,
                                   Origin,
                                   Species,
                                   Color) %>%
                     dplyr::distinct(),
                   by = c("label" = "Accession ID")) %>%
  dplyr::mutate(`Accession` = stringr::str_c(`label`,": ",Accession),
                  `Accession` = factor(`Accession`, levels = unique(`Accession`)))

dendro_plot <- ggplot() + 
  geom_segment(data = fit$segments,
               aes(x = x,
                   y = y,
                   xend = xend,
                   yend = yend)) +
  coord_flip() +
  scale_y_reverse(limits = c(1.0,-0.85),
                  breaks = seq(0,1,0.2)) +
  ylab("Scaled Distances")+
  ggplot2::theme(panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        axis.title.x = element_text(colour = NA), 
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
        )

The plots below each visualize the data clustered by \(\lvertΔm\rvert\)—the absolute difference in slopes—with label colors indicating either species, color, or origin. The distance measure (\(\lvertΔm\rvert\)) is standardized to be between zero and one.

As stated above, clustering by slope reveals a close relationship between the Dyck amaranth seeds, accession A (C02-007 Hopi Red Dye (C002)), and accession E (C02-013 Alegria (C008)). The former is intriguing in that it is geographically proximate to the Dyck cliff dwelling; the latter is interesting in that it shares its origin (Morelos, Mexico) with accession U, which clustered with the Dyck samples in the analysis above. Otherwise, there seems to be a slight propensity for accessions to cluster with those of similar color. The relationship between place of origin and slope cluster is inconclusive without further geographic analysis.

dendro_plot +
  geom_text(data = fit$labels,
            aes(x = x,
                y = y - 0.025,
                label = Accession,
                hjust = "outward",
                color = Species),
            size = 3) 

dendro_plot +
  geom_text(data = fit$labels,
            aes(x = x,
                y = y - 0.025,
                label = Accession,
                hjust = "outward",
                color = Color),
            size = 3) 

dendro_plot +
  geom_text(data = fit$labels,
            aes(x = x,
                y = y - 0.025,
                label = Accession,
                hjust = "outward",
                color = Origin),
            size = 3)